%func:@使用解析模型计算球形磁铁的磁感应强度(B)
%inpt:@M:剩磁场
%     @(ao,l):磁铁参数
%     @(xd,yd,zd)被计算的位置
%oupt:@B = [Bx;By;Bz]
%note:
%date:2021-07-04 by xiyin,li in Hust;
%@Copyright 2021 Xiyin,li.
function B = PMBall3D(M,R,position)
xd = position(1);
yd = position(2);
zd = position(3);
nu0 = 4*pi*10^(-7);
%% Method1:使用不推荐的积分函数，其更稳定
%% Method2:使用被推荐的积分函数，其速度更快
if(0)
    Bx1 = dblquad(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,0),0,R,0,(2*pi),1e-4)/(4*pi)*M;
    Bx2 = dblquad(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,1),0,R,0,(2*pi),1e-4)/(4*pi)*M;
    By1 = dblquad(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,2),0,R,0,(2*pi),1e-4)/(4*pi)*M;
    By2 = dblquad(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,3),0,R,0,(2*pi),1e-4)/(4*pi)*M;
    Bz1 = dblquad(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,4),0,R,0,(2*pi),1e-4)/(4*pi)*M;
    Bz2 = dblquad(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,5),0,R,0,(2*pi),1e-4)/(4*pi)*M;
else if(1)
        Bx1 = integral2(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,0),0,R,0,(2*pi),'RelTol',1e-4)/(4*pi)*M;
        Bx2 = integral2(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,1),0,R,0,(2*pi),'RelTol',1e-4)/(4*pi)*M;
        By1 = integral2(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,2),0,R,0,(2*pi),'RelTol',1e-4)/(4*pi)*M;
        By2 = integral2(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,3),0,R,0,(2*pi),'RelTol',1e-4)/(4*pi)*M;
        Bz1 = integral2(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,4),0,R,0,(2*pi),'RelTol',1e-4)/(4*pi)*M;
        Bz2 = integral2(@(r,theta)PMBallFunc(R,r,theta,xd,yd,zd,5),0,R,0,(2*pi),'RelTol',1e-4)/(4*pi)*M;
    end
end
B = [Bx1-Bx2;By1-By2;Bz1-Bz2];

%func:@使用解析模型计算球形磁铁时的被积函数
%inpt:@R:球的半径
%     @(r,theta):被积平面的参数方程
%     @(xd,yd,zd)被计算的位置
%     @flag |0   |1   |2   |3   |4   |5
%     @axis |x↑ |x↓ |y↑ |y↓ |z↑ |z↓
%oupt:@difB = dif[Bx;By;Bz]
%note:
%date:2021-07-04 by xiyin,li in Hust;
%@Copyright 2021 Xiyin,li.
function dif_B = PMBallFunc(R,r,theta,xd,yd,zd,flag)
%% 判断输入合法性
if(r>R)
    warnming("the r should NOT Greater than R")
end
%% 选择平面(上球面和下球面)
switch (mod(flag,2))
    case 0
        z = sqrt(R.^2-r.^2);
    case 1
        z = -sqrt(R.^2-r.^2);
end
x = r.*cos(theta);
y = r.*sin(theta);
R1 = sqrt((xd-x).^2+(yd-y).^2+(zd-z).^2);
%% 选择轴
switch (floor(flag/2)+1)
    case 1
        dif_B = r.*(xd-x)./R1.^3;
    case 2
        dif_B = r.*(yd-y)./R1.^3;
    case 3
        dif_B = r.*(zd-z)./R1.^3;
    otherwise
        warning("The flag is NOT allowed!");
end
